Spin-density-wave instability in graphene doped near the van Hove singularity 
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We study the instability of the metalUc state towards the formation of a new ground state in 
graphene doped near the van Hove singularity. The system is described by the Hubbard model and 
a field theoretical approach is used to calculate the charge and spin susceptibility. We find that 
for repulsive interactions, within the random phase approximation, there is a competition between 
ferromagnetism and spin-density wave (SDW). It turns out that a SDW with a triangular geometry 
is more favorable when the Hubbard parameter is above the critical value Uc{T), which depends on 
the temperature T, even if there are small variations in the doping. Our results can be verified by 
ARPES or neutron scattering experiments in highly doped graphene. 
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I. INTRODUCTION 

Graphene, a newly realized two-dimensional crystal of 
carbon atoms ordered on a honeycomb lattice^ is be- 
ing extensively studied due to its unusual electronic and 
structural properties. Undoped graphene is a zero-gap 
semiconductor with a linear low-energy dispersion rela- 
tion and a vanishing density of states at the Fermi level.^ 
Although the Coulomb interaction is unscreened in this 
regime, the system behaves mostly as a non-interacting 
electron liquid, where the minor effects of interactions 
are encoded in a renormalization of the Fermi velocity 
and the quasiparticle lifetime.^ The opposite case is that 
of a divergent density of states, the so called Van Hove 
(VH) singularity, which is associated to a saddle point in 
the band dispersion.l^ If the Fermi level lies near such a 
singularity, screening is perfect for wave- vectors connect- 
ing VH singularities and correlation effects may be en- 
hanced. Recently, several experimental groups have suc- 
ceeded on making this regime accessible for graphene .'^"^ 
The techniques used for this aim are different: in Ref. 5 
they have used twisted graphene, obtained from a rota- 
tion between stacked graphene layers, which allows one 
to tune the position of the VH singularity. Other meth- 
ods involve chemical doping of graphene^ or intercala- 
tion of gold clusters between the graphene layers.^ This 
opens up the fascinating possibility to study correlated 
electronic phases in this material, such as superconduc- 
tivity, charge-density wave (CDW) or spin-density wave 
(SDW). 

The proximity of the Fermi level to a VH singular- 
ity in most of the cuprate superconductors has trig- 
gered large efforts to understand the role of a peaked 
density of states (DOS) on the electronic properties 
of an electron liquid.^ Superconductivity, itinerant fer- 
romagnetism (FM), CDW, and SDW are examples of 
competing instabilities associated with the VH sce- 
nario. In graphene, the existence of some of these in- 
stabilities has already been evidenced by recent experi- 
ments, such as CDW^ ^ or superconducting pairing due 



to electron-electron interactions,^ following the Kohn- 
Luttinger mechanism.^ CDW and SDW phases usu- 
ally occur in systems with Fermi surface nesting, i.e. 
when the Fermi surface can be mapped onto itself by 
a (nonzero) k- vector (nesting vector). At low doping, 
graphene shows some nesting, since the Fermi surface 
is composed of two circles around the K and K' points. 
However, since the DOS is very low, these peaks are small 
and a very high (Hubbard) interaction is required to enter 
the DW regime.!^ The situation changes if we strongly 
dope graphene around the value /i ^ where fi is the 
chemical potential and t ^ 2.8eV is the hopping param- 
eter. The Fermi surface then acquires a hexagonal shape 
and nesting occurs in three directions (see Fig. [T]). Fur- 
thermore, the VH singularities lie at the Fermi surface 
for this doping. At these points, the DOS diverges and 
we expect some non-trivial peaks in the susceptibility. 

In this paper, we investigate the Hubbard model in a 
honeycomb lattice using a field theoretical approach. By 
performing a Hubbard- Stratonovich transformation, we 
determine the effective action in terms of an 8-component 
order parameter, which accounts for charge- and spin- 
degrees of freedom in the A and B graphene sublattices. 
We find that the charge susceptibility never diverges for 
repulsive interactions, thus excluding the possibility of a 
CDW formation, whereas the spin susceptibility exhibits 
several peaks. The peaks at the nesting wave vectors k 
are stronger than the one at k = (see Fig. 2|, thus 
signaling that SDW is the leading instability, that wins 
against FM for repulsive interactions. 

In the following, we introduce the model and outline 
the main steps of the calculations in Sec. II. Then, we 
present the temperature vs interaction vs doping phase 
diagram and show that the metallic state becomes un- 
stable towards a SDW phase for realistic values of the 
interaction parameter U in Sec. HI. Finally, we discuss 
possible experimental techniques which could be used to 
observe our results and draw our conclusions in Sec. IV. 
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Figure 1: (Color online) (a): Blue dashed line is the Brillouin 
zone of monolayer graphene. Red line is the Fermi surface for 
\i — t. Red dots are the van Hove singularities and the green 
arrows are the nesting vectors. 



where the action is given by 
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In this expression, r is the imaginary time variable, {3 = 
l/ksT and H is defined in Eq. ([T]). We will investigate 
the possible appearance of CDW and SDW instabilities 
by using a path integral formalism recently developed by 
some of the authors.^ In the following, we briefly outline 
the procedure. 

We start by performing a Hubbard- Stratonovich trans- 
formation that eliminates the quartic term in the action, 
but introduces eight auxiliary bosonic fields p^, p^, M^, 
and related, respectively, to the electronic and the 
magnetic degrees of freedom of the fermionic fields of 
each sublattice. In Fourier space, the action then reads 



II. THE MODEL 

Due to the strong screening of interactions by the elec- 
tron gas at the VH filling, we consider only the local 
Coulomb repulsion (the usual Hubbard U term), which 
is a good approximation around this doping level.l^^inSI 
This term leads to a deformation of the band disper- 
sion toward the saddle point an e^ffect that has been 
observed experimentally by ARPES.^ Notice that longer- 
range interaction terms, such as the nearest-neighbor re- 
pulsion, would lead to a richer phase diagram, with the 
possibility of a Pomeranchuk instability.^^ This sce- 
nario is in sharp contrast with low doped graphene, for 
which the DOS vanishes and the Coulomb interaction 
will be only slightly screened and therefore will be long 
ranged. Hence, we use a tight-binding model with a Hub- 
bard interaction /7, 

H = -t J2 (4,shs + H.c.) + U J2 4,tS\4^i.^^i.t, 

<i,j>,s jeA,B 

(1) 

where the operator Cj can be either aj or 6j, depending 
whether j is a label of the A or B sublattice, respec- 
tively. Defining Cj = (cj^,c^^), the Hubbard term can 
be rewritten using the relation. 
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S • S , 



where rij = CjCj is the on site number operator and Sj = 

(1/2) Cj a Cj, describes the spin on the lattice site j. Note 
that the inner product Sj - A, where A is some vectorial 
field, is defined by Sj ■ A = (1/2) Cj (cr • A) Cj, and that 
a is the vector of Pauli matrices. The grand-canonical 
partition function describing the system reads 



J d[c^]d[c]e" 



-S[c\c]/h 



5 = -^ E (<ni„)-Gk,U',n' 



k,n,k',n' 



where a = a,b. The inverse Green's function is defined 



^k,n;k',n' ~ ^Ok,n;k',n' 

Green's function reads 



"^^Ok,n;k',n' 



-{jii^ihujn)! -^Tkl 



with I a 2 X 2 identity matrix, uj^ = 7r(2n + 1)/^/^ 
the fermionic Matsubara frequency, and 7k = e*^°^^ + 
e~*^°^y/^ cos(V3aoA:a,/2). In what follows we set the lat- 
tice constant ao = 1. The self energy is given by 



-1 



a-M^ - 



where N is the number of sites of each sublattice and 
both, and carry a subindex (k — k',n — n'). By 
introducing the eight component vector 

and integrating out the fermionic fields, we obtain the 
effective action 



77 • M_ 



/iTr[ln(-G-^)], 



where rj = diag(— 1, 1, 1, 1, — 1, 1, 1, 1) and the partition 
function becomes Z = j d[M] exp(— S'eff/fi). Finding the 
susceptibility from here is a two steps process. First, 
we introduce external source fields J and then we ex- 
pand the fields in the action around their mean field 
value Mk^^ = (Mk,n)j + ^^Mk,n- This results in a self- 
consistent equation for the mean field values of the fields. 
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Figure 2: (Color online) Density plot of the largest eigenvalue 
of TJX for T — O.Olt. Brighter regions correspond to higher 
peaks in the spin susceptibility. The peaks at the three nesting 
wavevectors (see Fig. [T]) are higher than the one at k = 0, 
indicating that the SDW instability wins over the FM one. 



Using this equation, we may determine how the bosonic 
fields react on a distortion of the source fields. 

The mean field values are those at which the action 
has an extremum. Thus, the linear terms in (^Mk^^ are 
imposed to vanish, yielding the self-consistent equation 



U 



p,m;p— k,m— nj? 



where P = [diag(I, 0), diag((j^, 0), diag((j^, 0), diag(cr^, 0), 
diag(0, 1), diag(0, cr^), diag(0, cr^), diag(0, cr^)]^ and the 
trace is taken over the matrices only, because the sum 
takes care of the trace in k- and frequency-space. After 
differentiating this equation with respect to the source 
J, one finds the susceptibility through the relation 



d(M)j = hUxT^dJ, 



(2) 



where 



Moreover, we can also determine the free energy by 
performing a Legendre transformation /3F[(M)j] = 
[/~^(M) J • J — ln(Z[J]) on the partition function which, 
up to quadratic order in the deviation A(M) j = (M) j — 
(M)o and without an additive constant, is 



2hm 



(A{M)^)t.(Xo^P^)-i.A{M)^. (3) 



The susceptibility Xo^^^ is evaluated in the absence of 
the source field, J = 0. From here on, we assume 
that the system is in a homogeneous state, for which 
the mean field values of the boson fields are given by 




Figure 3: (Color online) Plot of the largest eigenvalue of r]x 
for T = 0.0025t for kx = 0. Notice that the peak correspond- 
ing to SDW at ky = 27V is larger than the FM peaks at /c^ = 
and ky = 4:7r. 



(M^^)o = (Mq o)(^k,o^^n,o- As a consequence, the self- 
energy, the Green's function and the susceptibility matrix 
Xo are all diagonal in momentum and frequency space. 
Due to a nonzero self-energy, the Hamiltonian gets renor- 
malized, such that 



hG- 



(/i + ihujri)'i- - Ho 



(/i + ihWn)! — i^r( 



where Hq is given by the first term in Eq. (1). Assuming a 
(renormalized) Hamiltonian, which may be diagonalized 
by using the unitary operators U\^, such that Ul^H\JJ\^ = 
where /^^^ = dmg{Sa,i, Sa,2, Sa,3, Sa,^) and 



11^ the susceptibility can be written a; 



^P+k,p 



N 



tp 



rpT^r ;Q!,/3 
P+k,p ' 



E 



where Vt^ are bosonic Matsubara frequencies and 
npix) = l/(e^^ + 1) is the Fermi distribution function. 
Notice that the overlap between the electron and hole 
wave-functions is taken into account through the term 

rTir,r' •,a,(3 
^P+k,p • 



T/t(10-^) 


T(K) 


ur 




0.025 


0.87 


1.0570 


0.8444 


0.05 


1.74 


1.6710 


0.9954 


0.075 


2.6 


2.0410 


1.0776 


0.15 


5.2 


2.5018 


1.2184 


0.25 


8.7 


2.6828 


1.3322 


0.5 


17 


2.8896 


1.5118 


1 


35 


3.1272 


1.7288 


2.5 


87 


3.5082 


2.0920 


5 


175 


3.8652 


2.4444 


10 


350 


4.3030 


2.8878 



Table I: Numerical values of the critical coupling Uc (in units 
oft) for a FM and a SDW phases calculated on a 10000 x 10000 
mesh. 
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Figure 4: (Color online) (a) Plot of the critical coupling required to enter the SDW regime (red line) and the FM regime (blue 
line) as a function of temperature, starting from the metallic phase at small U. We observe that the SDW phase wins for all 
temperatures, (b) Zoom-in of figure (a) for low temperatures. Both figures have been determined using a 10000x10000 mesh, 
for which convergence was eventually reached for T > 0.5 x 10~^. (c) Same as (a), but calculated in a 700x700 mesh. When 
not enough points are taken into account, the two lines cross, thus leading to an erroneous conclusion that a FM phase would 
be more favorable at low temperatures. 



III. NUMERICAL RESULTS 

Next, we investigate the possible instabilities by con- 
sidering the static susceptibility %k(0)- To evaluate the 
latter, one needs to perform a sum over all momenta in 
the Brillouin zone. We proceed by dividing the Brillouin 
zone in an X mesh, which implies that the sum in the 
expression for the susceptibility has N'^ terms. The in- 
stability condition for repulsive interactions requires the 
interaction strength U to exceed the critical value Uc de- 
fined by = det {x J =Q — Ucr]) = det{r]XjLo~^c)' This re- 
lation links the critical interaction strength to the largest 
eigenvalue Ak of the matrix ?7%k(0) by 



u: 



maxAk = A 

k 



Q- 



In Fig. [2] a density plot of the largest eigenvalue of 
the matrix r7Xk(0) is shown for a temperature of T = 
0.0 It ^ 350 K. One can distinguish four inequivalent 
peaks (bright regions in Fig. [2]). The one at k = 
corresponds to a FM instability, while the other three, 
at nonzero k, correspond to the three independent nest- 
ing vectors (see Fig. [T]) and hence give rise to a SDW. 
Although not expected, for all temperatures the SDW 
peaks turn out to be higher than the peak correspond- 
ing to the FM phase (see Fig. 3|. However, this result 
depends heavily on the size of the mesh: for low temper- 
atures, the sum converges slowly, and a wrong result can 
be inferred if is not large enough. Indeed, by consider- 
ing a 700 x 700 mesh, we find t hat a FM phase would set 
in at low temperatures [see Fig. |4(c)] , whereas the results 
for a finer mesh (10000 x IQQQQ ) ind icate th at the true 
ground state is a SDW [see Fig. |4(a)| and Fig. |4(b)| . 

For zero temperature and a doping value exactly at 
the VH singularity, the critical coupling for both FM and 
SDW is zero. As explained below, our formalism is not 
suited to determine which state is more favorable. For fi- 
nite temperatures, the SDW phase transition has a lower 
critical U than the FM one. Therefore, between the red 
and blue lines in Fig. [4] we find a SDW phase. Above the 



blue line, the system allows for a FM phase transition. 
The starting point of our formalism is a homogeneous 
ground state, after which we can determine the critical 
U that is needed to enter a more ordered phase, like the 
FM, SDW or CDW phase. To reach the regime above the 
blue line in Fig. |4j one has to start from a SDW ground 
state and subsequently determine if a phase transition to 
the FM phase will occur. Since the SDW ground state 
is inhomogeneous, our formalism is thus unable to deter- 
mine the leading instability when U is in this regime. 

If the system is tuned away from optimal doping, the 
VH singularities will no longer lie on the Fermi surface 
and nesting is reduced in general. This will result in 
a lowering of the peaks in the maximal eigenvalue of 
the susceptibility and hence the critical couplings will 
increase. This behavior is shown in Fig. [5j The SDW 
peak is always higher than the FM one, but the height 
decreases quite rapidly as function of doping. The inter- 




Figure 5: Largest eigenvalue vs temperature and vs chemical 
potential. The critical U scales as 1/Aq and will therefore 
increase if the doping level is tuned away from /x = \t\. 
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action strength in graphene cannot be tuned externally. 
Although at low doping the onsite Hubb ard interaction 
is estimated to be around U/t = 3.5121111 values of U for 
highly doped graphene are, to the best of our knowledge, 
not yet known, but should not be different from the low- 
doping values. We therefore expect that a SDW phase 
should be observed experimentally, in the regime of T 
and U parameters shown in Fig. [s] (see also Table 

IV. DISCUSSIONS AND CONCLUSIONS 

Recently, a CDW has been experimentally observed 
in twisted graphene bilayers doped at the VH singular- 
ity, but it has disappeared when the hopping parameter 
between the two layers has been put to zero.^ A CDW 
phase seems also to be present for graphene grown on 
top of a superlattice of gold intercalated clusters,^ and a 
superconducting phase has been conjectured to occur for 
graphene heavily doped with Ca and K on both sides. ^ 
However, it is unclear which is the appropriate theoretical 
model to describe graphene under these circumstances, 
and for the conventional Hubbard model in a single layer, 
a CDW can only occur for attractive on-site interactions. 

In conclusion, we have evaluated possible SDW and 
CDW instabilities from a metallic phase for the honey- 
comb lattice of a single layer of graphene doped at the 



VH filling. We found that charge and spin degrees of free- 
dom are decoupled, and that CDWs occur for attractive, 
whereas SDWs occur for repulsive interactions. A peak in 
the spin-susceptibility at zero wave- vector has also been 
found, although precise calculations indicate that it is al- 
ways sub-dominant with respect to the ones appearing at 
the nesting wavevectors, thus determining that the SDW 
instability must win over FM in the neighborhood of the 
metallic phase. This finding contradicts a previous re- 
sult in the literature, and could only be reached after a 
very careful and time consuming computation in a very 
large mesh. An interesting point is that the SDW phase 
is expected to occur at rather high temperatures. This is 
not very surprising, given that phenomena such as CDW 
or the integer quantum Hall effect, that usually occur 
at very low temperatures, have already been detected at 
high temperatures in graphene (T = 77K for the first and 
room temperature for the latter).^ We hope that our 
results may stimulate further experiments to unveil the 
existence of a SDW phase in heavily doped graphene. 
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